#version 120
#extension GL_EXT_geometry_shader4 : enable
#extension GL_EXT_gpu_shader4 : enable

uniform mat4 Q0;
uniform mat4 Q1;
uniform mat4 Q2;
uniform mat4 Q3;

mat4 Q[4];
varying in vec4 osTetra[4];

varying out mat4 Qi;
varying out mat4 Qo;
varying out vec4 pi;
varying out vec4 po;
varying out vec4 cori;
varying out vec4 coro;

ivec4 sort(void)
{
	vec3 n = vec3(1.0,0.0,0.0);
	ivec4 idx = ivec4(0,1,2,3);
	vec4 coss;
	int t;
  
  //<MaisAbaixo>
	t = 0;
	for(int i = 1; i < 4; ++i)
		if(gl_PositionIn[i].y < gl_PositionIn[t].y)
			t = i;
  
	idx[0] = t;
	idx[t] = 0;
  //</MaisAbaixo>
	
	coss[0] = dot(n, normalize(vec3(gl_PositionIn[0].xy,0.0) - vec3(gl_PositionIn[idx[0]].xy,0.0)));
	coss[1] = dot(n, normalize(vec3(gl_PositionIn[1].xy,0.0) - vec3(gl_PositionIn[idx[0]].xy,0.0)));
	coss[2] = dot(n, normalize(vec3(gl_PositionIn[2].xy,0.0) - vec3(gl_PositionIn[idx[0]].xy,0.0)));
	coss[3] = dot(n, normalize(vec3(gl_PositionIn[3].xy,0.0) - vec3(gl_PositionIn[idx[0]].xy,0.0)));
  
	if(coss[idx[2]] > coss[idx[1]] && coss[idx[2]] > coss[idx[3]])
	{
		t = idx[1];
		idx[1] = idx[2];
		idx[2] = t;
	}else if(coss[idx[3]] > coss[idx[1]] && coss[idx[3]] > coss[idx[2]])
	{
		t = idx[1];
		idx[1] = idx[3];
		idx[3] = t;
	} // else idx[1] ja eh o menor
  
	if(coss[idx[3]] > coss[idx[2]])
	{
		t = idx[2];
		idx[2] = idx[3];
		idx[3] = t;
	} // else idx[2] ja eh o menor
  
	return idx;
  
}

mat3 inverse3( mat3 A )
{
  mat3 invA ;
  float d = +A[0][0]*(A[1][1]*A[2][2]-A[2][1]*A[1][2])
            -A[0][1]*(A[1][0]*A[2][2]-A[1][2]*A[2][0])
            +A[0][2]*(A[1][0]*A[2][1]-A[1][1]*A[2][0]);
  float invdet = 1.0/d;
  invA[0][0] =  (A[1][1]*A[2][2]-A[2][1]*A[1][2])*invdet;
  invA[1][0] = -(A[0][1]*A[2][2]-A[0][2]*A[2][1])*invdet;
  invA[2][0] =  (A[0][1]*A[1][2]-A[0][2]*A[1][1])*invdet;
  invA[0][1] = -(A[1][0]*A[2][2]-A[1][2]*A[2][0])*invdet;
  invA[1][1] =  (A[0][0]*A[2][2]-A[0][2]*A[2][0])*invdet;
  invA[2][1] = -(A[0][0]*A[1][2]-A[1][0]*A[0][2])*invdet;
  invA[0][2] =  (A[1][0]*A[2][1]-A[2][0]*A[1][1])*invdet;
  invA[1][2] = -(A[0][0]*A[2][1]-A[2][0]*A[0][1])*invdet;
  invA[2][2] =  (A[0][0]*A[1][1]-A[1][0]*A[0][1])*invdet;
  
  return invA ;
}

/*vec3 parametroInterpolacao2(in ivec3 idx,in vec2 p)
{
	mat3 A = mat3(vec3(gl_PositionIn[idx[0]].xy, 1.0),  //coluna
                vec3(gl_PositionIn[idx[1]].xy, 1.0),  //coluna
                vec3(gl_PositionIn[idx[2]].xy, 1.0)); //coluna
	vec3 P = vec3(p,1.0);
  
	return inverse3(A)*P;
}*/
 // Quero verificar se o algorítimo de baixo e mais rapido que o de cima
 vec3 parametroInterpolacao2(in ivec3 idx,in vec2 p)
 {
    float areaTotal;
    vec3 area;
    vec3 p0 = vec3(gl_PositionIn[idx[0]].xy, 0.0);
    vec3 p1 = vec3(gl_PositionIn[idx[1]].xy, 0.0);
    vec3 p2 = vec3(gl_PositionIn[idx[2]].xy, 0.0);
    vec3 pp = vec3(p, 0.0);

    areaTotal = length(cross(p1-p0, p2-p0));// /2.0;

    area[0] = length(cross(p2-p1,pp-p1));// /2.0;
    area[1] = length(cross(p0-p2,pp-p2));// /2.0;
    area[2] = length(cross(p1-p0,pp-p0));// /2.0;

    area[0]/=areaTotal;
    area[1]/=areaTotal;
    area[2]/=areaTotal;

    return area;
 }

vec2 intersection(in ivec4 idx)
{
	vec2 ret = vec2(0.0);
	vec2 a = gl_PositionIn[idx[0]].xy;
	vec2 b = gl_PositionIn[idx[2]].xy;
	vec2 c = gl_PositionIn[idx[1]].xy;
	vec2 d = gl_PositionIn[idx[3]].xy;
  
	float der = (b.y - a.y)/(b.x - a.x);
	float ydc = d.y-c.y;
	float xdc = d.x-c.x;
  
	ret.t = (a.y - c.y + der*(c.x-a.x)) / (ydc-der*xdc);
	ret.s = (c.x + ret.t*xdc - a.x)/ (b.x - a.x);
  
	return ret;
}


void main()
{
	ivec4 idx;
	vec2 s,t;	
	vec4 center;
	mat4 Qic;
	mat4 Qoc;
	vec4 pic;
	vec4 poc;
	vec4 coric;
	vec4 coroc;
	int nTriangulos;

	Q[0] = Q0;
	Q[1] = Q1;
	Q[2] = Q2;
	Q[3] = Q3;  

	idx = sort();
	
	vec2 p = intersection(idx);
	s[0] = 1.0 - p.s;
	s[1] = p.s;
	t[0] = 1.0 - p.t;
	t[1] = p.t;
  
	if(s[0] > 0.0 && s[0] < 1.0 && t[0] < 1.0 && t[0] > 0.0)
	{
		float z0 = s[0]*gl_PositionIn[idx[0]].z + s[1]*gl_PositionIn[idx[2]].z;
		float z1 = t[0]*gl_PositionIn[idx[1]].z + t[1]*gl_PositionIn[idx[3]].z;
		if(z0 < z1)
		{
			center = s[0]*gl_PositionIn[idx[0]] + s[1]*gl_PositionIn[idx[2]];
			
			pic =    s[0]*osTetra[idx[0]] + s[1]*osTetra[idx[2]];
			Qic =    s[0]*Q[idx[0]]          + s[1]*Q[idx[2]];
			coric =  s[0]*gl_FrontColorIn[idx[0]]     + s[1]*gl_FrontColorIn[idx[2]];
      
			poc =    t[0]*osTetra[idx[1]] + t[1]*osTetra[idx[3]];
			Qoc =    t[0]*Q[idx[1]]          + t[1]*Q[idx[3]];
			coroc =  t[0]*gl_FrontColorIn[idx[1]]     + t[1]*gl_FrontColorIn[idx[3]];
			
		}else
		{
			center = t[0]*gl_PositionIn[idx[1]] + t[1]*gl_PositionIn[idx[3]];
			
			pic =    t[0]*osTetra[idx[1]] + t[1]*osTetra[idx[3]];
			Qic =    t[0]*Q[idx[1]]          + t[1]*Q[idx[3]];
			coric =  t[0]*gl_FrontColorIn[idx[1]]     + t[1]*gl_FrontColorIn[idx[3]];
			
			poc =    s[0]*osTetra[idx[0]] + s[1]*osTetra[idx[2]];
			Qoc =    s[0]*Q[idx[0]]          + s[1]*Q[idx[2]];
			coroc =  s[0]*gl_FrontColorIn[idx[0]]     + t[1]*gl_FrontColorIn[idx[2]];
		}
		nTriangulos = 4;
	}else
	{
		int t = idx[2];
		idx[2] = idx[3];
		idx[3] = t;
		// idx[0,1,2] formam um triangulo CCW
    
		vec3 pI = parametroInterpolacao2(idx.xyz, gl_PositionIn[idx[3]].xy);
    
		float z0 = pI[0]*gl_PositionIn[idx[0]].z + pI[1]*gl_PositionIn[idx[1]].z + pI[2]*gl_PositionIn[idx[2]].z; //z do triangulo
		float z1 = gl_PositionIn[idx[3]].z;
		if(z0 < z1) // triangulo na frente
		{
			center = pI[0]*gl_PositionIn[idx[0]] + pI[1]*gl_PositionIn[idx[1]] + pI[2]*gl_PositionIn[idx[2]]; 
			
			pic =    pI[0]*osTetra[idx[0]] + pI[1]*osTetra[idx[1]] + pI[2]*osTetra[idx[2]];
			Qic =    pI[0]*Q[idx[0]]          + pI[1]*Q[idx[1]]          + pI[2]*Q[idx[2]];
			coric =  pI[0]*gl_FrontColorIn[idx[0]]     + pI[1]*gl_FrontColorIn[idx[1]]     + pI[2]*gl_FrontColorIn[idx[2]];
      
			poc =    osTetra[idx[3]];
			Qoc =    Q[idx[3]];
			coroc =  gl_FrontColorIn[idx[3]];
		}else // vertice na frente
		{
			center = gl_PositionIn[idx[3]];
      
			pic =    osTetra[idx[3]];
			Qic =    Q[idx[3]];
			coric =  gl_FrontColorIn[idx[3]];
      
			poc =    pI[0]*osTetra[idx[0]] + pI[1]*osTetra[idx[1]] + pI[2]*osTetra[idx[2]];
			Qoc =    pI[0]*Q[idx[0]]          + pI[1]*Q[idx[1]]          + pI[2]*Q[idx[2]];
			coroc =  pI[0]*gl_FrontColorIn[idx[0]]     + pI[1]*gl_FrontColorIn[idx[1]]     + pI[2]*gl_FrontColorIn[idx[2]];
		}
    
		nTriangulos = 3;
	}
  
	gl_FrontColor = vec4(1.0);
	for(int i = 0; i < nTriangulos; i++)
	{
		int iplus1 = (i+1)%nTriangulos;
		gl_Position = center;
		Qi = Qic;
		Qo = Qoc;
		pi = pic;
		po = poc;
		cori = coric;
		coro = coroc;
		EmitVertex();
    
		gl_Position = gl_PositionIn[idx[i]];
		Qi = Q[idx[i]];
		Qo = Q[idx[i]];
		pi = osTetra[idx[i]];
		po = osTetra[idx[i]];
		cori = gl_FrontColorIn[idx[i]];
		coro = gl_FrontColorIn[idx[i]];
		EmitVertex();
    
		gl_Position = gl_PositionIn[idx[iplus1]]; 
		Qi = Q[idx[iplus1]];
		Qo = Q[idx[iplus1]];
		pi = osTetra[idx[iplus1]];
		po = osTetra[idx[iplus1]];
		cori = gl_FrontColorIn[idx[iplus1]];
		coro = gl_FrontColorIn[idx[iplus1]];
		EmitVertex();
    
		EndPrimitive();
	}
}
